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Abstract 



We report the results of simulations of the Lebwohl-Lasher model of the 
nematic-isotropic transition using a new cluster Monte Carlo algorithm. The 
algorithm is a modification of the Wolff algorithm for spin systems, and greatly 
reduces critical slowing down. We calculate the free energy in the neighbor- 
hood of the transition for systems up to linear size 70. We find a double well 
structure with a barrier that grows with increasing system size, obeying finite 
size scaling for systems of size greater than 35. We thus obtain an estimate 
of the value of the transition temperature in the thermodynamic limit. 
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The Lebwohl-Lasher (LL) model JT[] is a lattice model of rotors with an orientational 
order-disorder transition. The long axes of the rotors are specified by unit vectors <jj. 
While it neglects the coupling between the orientational and translational degrees of freedom 
present in a real nematic liquid crystal, it is generally believed that this coupling does 
not play a significant role at the nematic-isotropic (NI) transition. With the absence of 
translational degrees of freedom the LL model is particularly well-suited for large-scale 
simulations of the transition. The model is defined by the Hamiltonian: 



where the sum is over all nearest-neighbors and e is a coupling parameter. The LL model 
has been intensively investigated using Monte Carlo techniques since its introduction 0-0. 

As in real experimental systems the NI transition in the LL model is weakly first-order; 
thus, there is significant critical slowing down in the neighborhood of the transition. In a 
Monte Carlo simulation the system gets trapped in one of the free energy wells corresponding 
to the nematic or isotropic phase, and the conventional single flip Metropolis algorithm be- 
comes inefficient especially as the system size is increased. While Boschi et al. |7| carried out 
simulations on systems as large as 120 x 120 x 120, the most detailed study of the NI transi- 
tion was carried out by Zhang et al. [g] on systems up to 28 x 28 x 28. These authors used the 
Lee-Kosterlitz finite size scaling method |||§, supplemented by the Ferrenberg-Swendsen 
reweighting technique 0] to determine the order of the NI transitions and estimate the 
value of the transition temperature T c in the thermodynamic limit. In the Lee-Kosterlitz 
method one examines the finite size scaling of the free energy barrier AF between the ne- 
matic and isotropic phases; at a first-order transition this should be an increasing function 
the linear system size L, while it should approach a constant for systems with continuous 
phase transitions. For a large enough system, specifically L ^> £, where £ is the correlation 
length, a finite size scaling analysis predicts that AF ~ L 2 for three-dimensional systems. 
In the LL model Zhang et al. found a small free energy barrier appearing at the two largest 
system sizes they studied, L = 24 and 28, and thus did not have enough data to carry out 




(1) 
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a finite size scaling analysis of AF. Instead they estimated the value of T c in the thermo- 
dynamic limit by extrapolating three different measures of T c : the positions of the maxima 
in the specific heat and susceptibility and the the temperature where the two free energy 
wells are of equal depth. However, as we demonstrate below, the system sizes considered by 
Zhang et al. are not in the finite-size scaling regime, and thus their estimate of T c in the 
thermodynamic limit is not accurate. 

Over the past decade significant advances have been made in algorithm development 
which overcome critical slowing down in magnetic spin systems [11]. In particular single 



cluster algorithms have proven to be very efficient in simulating the three-dimensional Ising, 
XY and Heisenberg models. These algorithms are nonlocal updating methods where a single 
cluster of spins is constructed and the spins within the cluster are updated simultaneously. 
In the Ising case |12| clusters of spins are formed by creating bonds between parallel spins 



with a probability that guarantees detailed balance. For models with continuous symmetry 
Wolff jnj introduced a cluster algorithm where "parallel spins" refer to spins which point 
in the same hemisphere. A hemisphere is defined by an equatorial plane perpendicular 
to a randomly chosen direction f . Nematic liquid crystals differ in an important symmetry 
aspect from magnetic systems, namely, "up" and "down" spins are equivalent. To construct a 
cluster algorithm suitable for simulating the LL model we have modified the Wolff algorithm 
to account for this symmetry difference. As in the original Wolff algorithm we randomly 
choose a direction f. Then we reflect any molecular long axes for which oi ■ f < 0, by the 
transformation a — > —a; note that the Hamiltonian Ti is invariant under this operation. 
Next we choose a site i at random and reflect it, a[ = i2(?)oi using the reflection operator 
R(r) defined by: 

R(r)ai = -a { + 2{<n ■ r)r. (2) 

This operation is illustrated in Fig. |l]a. Unlike the original Wolff reflection operation which 
reflects spins from the original hemisphere to the opposite one, the present reflection operator 
keeps the molecular orientation vectors in the same hemisphere defined by f . Next we form 
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bonds with the nearest-neighbors of a\ with probability: 

Pii = 1 - expjmiJo,/^ • a 3 f - (3(a' t ■ (i?(r>,)) 2 



1 — exp jmin 



0,4/%; • v){a 3 ■ v)({a[ ■ a 3 ) - {a[ ■ t)(a s • ?))] | 



(3) 



where (3 = 3e/2ksT. This probability is a modification of the one introduced by Wolff, 
replacing the Heisenberg interaction — Jovcr,- by the Lebwohl-Lasher interaction — -e^-Cj) 2 . 
As in the original Wolff algorithm we continue this process, forming bonds with the nearest- 
neighbors of all reflected molecular orientation vectors until the cluster cannot grow any 
further. 

To understand the formation of clusters, consider the projection of two molecular orien- 
tation vectors, <7j and one of its nearest-neighbors <jj, on the plane perpendicular to f before 
the reflection operation is performed (see Fig. [l|b). A bond between these two molecules 
will likely form if the angle between their projections is less than 90°. Note that the 
probability Eq. (Q) for the bond formation is maximized when the angle between a\ and a 3 is 
90° and each of these molecules makes an angle of 45° with f . Thus, as in the original Wolff 
algorithm at low temperatures the molecules are nearly all aligned and it is highly probable 
that large fraction of all molecules will be flipped at once. On the other hand at high tem- 
peratures the distribution of molecules will be isotropic, resulting in flipping small clusters 
in random directions. In the intermediate region close to the NI transition temperature the 
system flips between isotropic and nematic states. 

We have used our cluster algorithm to simulate the LL model on a simple cubic lattice of 
linear dimension L, 30 < L < 70, with periodic boundary conditions, in order to study the 
properties of the NI transition. The temperature T was measured in dimensionless units of 
e/fce, in agreement with the units used in previous studies of this model. Our initial random 
configurations were equilibrated at least 200000 Monte Carlo steps (MCS) before starting 
production runs. We found that the average cluster size at temperatures close to the NI 
transition is approximately 0.17iV sites per cluster (where N is the total number of lattice 
sites), essentially independent of system size. Approximately half of the MCS resulted in 



5 



clusters with fewer than ten sites and were efficiently simulated with scalar code. However, a 
significant fraction of clusters had N/2 sites or greater, and employing a vectorizable cluster 
construction method |nj yielded a sixfold speedup. 

For each configuration generated, we calculated the energy per site, E — 71/ N. To 
ascertain the nature of the phase transition we proceeded as in Ref. || and used the method 
of Lee and Kosterlitz HH, which relies on the single histogram reweighting technique of 
Ferrenberg and Swendsen flO |. Following the approach of the latter authors we stored 



the configuration data in a histogram H(E,T,L). The normalized probability distribution 
function P(E, T, L) of the energy is then given by: 

Given this distribution function at temperature T, the Ferrenberg-Swendsen method allows 
the calculation of thermodynamic quantities at a different temperature T' in the neigh- 
borhood of T. Specifically, thermodynamic quantities at T' can be calculated using the 
distribution function P(E, T', L) where: 

, H(E^L)exp(-A(3E) 

P{E ' T ' L) -Y, E H(E,T,L)eM-WE) (5) 

and 

A(3= (l/T - 1/T). (6) 

Thus, accurate information over the entire critical region can be extracted from a small 
number of simulations. 

The Lee-Kosterlitz method utilizes the system size dependence of the barrier AF separat- 
ing the isotropic and nematic free energy minima at the transition temperature to determine 
the order of the transition. If the barrier grows with increasing L then the transition is first 
order; furthermore, if finite size scaling holds, then AF ~ L 2 in a three-dimensional system. 
To determine the barrier height we use the free-energy-like quantity : 

F(E,T,L)~-lnP(E,T,L) (7) 
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which differs from the true free energy by additive quantities dependent only on T and 
L which are irrelevant to computations of free energy differences. This free-energy-like 
quantity is shown in Fig|| for different system sizes, and we note the appearance of a 
pronounced double-well structure for sufficiently large system sizes. The right and left hand 
wells correspond to the nematic and isotropic phases respectively. In collecting our data we 
made sure that system made at least 100 hops between the two wells for the largest system 
size of 70, for a run of 6 x 10 6 MCS. For each system size, we performed sufficient MCS such 
that the typical number of points in one bin of the histogram H is much larger than the 
variation in the exponential factor exp(ATL 3 ), where AT — T' — T. This criterion arises 
from the requirement that the peak in the reweighted distribution Eq. (^) avoid the "wings" 
of the measured histogram. Typically we found approximately 10 4 points in each bin and 
the exponential factor varied by about 10. 

Zhang et al. § made similar plots of the free energy F(S, T, L) as a function of the 
nematic order parameter S rather than the energy E. For the system sizes studied by these 
authors, with L < 28, the free energy function F(E, T, L) is a much weaker indicator of the 
nature of the NI transition. However, for the system sizes we have studied, with 30 < L < 70, 
the free energy as a function of E is a very good indicator as illustrated in Fig.|2|. We 
calculated F(S,T,L) for L = 28, the largest system size studied by Zhang et al. to check 
that our cluster algorithm yields the same transition temperature they determined using 
the conventional single spin flip MC algorithm. In general we did not calculate F(S, T, L) 
because this requires calculation of a histogram H(E, S, T, L), dependent on 5* as well as E 
in order to carry out the reweighting. Calculation of this multiple histogram with sufficiently 
good statistics is prohibitively time consuming for large systems. 

The barrier height AF(L) can be computed as follows : 

AF(L) = F(E mj T, L) - F(E U T, L) (8) 

where E m is the energy corresponding to the top of the free energy barrier and E% is either 
one of the degenerate local minima. Our results for the barrier height as a function of system 
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size are shown in Fig.^, where we see that finite size scaling holds for systems of size L > 35, 
i.e. beyond the largest systems considered in ref. ||. The transition temperature T C (L) for a 
particular system size L is given by the value of the temperature where the two free energy 
wells have equal depths. Our results for T C (L) are shown in Fig.^. From the straight line 
plotted in the figure we see that the finite size scaling relation, 

AT = T C -T C (L) ~L-\ (9) 

works well for systems of size L ^ 35, as we would expect, given the scaling of the free 
energy barrier, AF(L). The intersection of the straight line with the T C (L) axis yields our 
estimate of the transition temperature in the thermodynamic limit: T c = 1.1225 ± 0.0001. 
This value is lower than that obtained in Ref. ||, (T c = 1.1232±0.0001). However, our data 
indicates that the system sizes studied in Ref. were not large enough to carry out a finite 
size scaling analysis. Even though the three measures of T C (L) used in this latter reference 
apparently extrapolate to a single number, there is no justification for using a straight line 
extrapolation for the transition temperature in the absence of finite size scaling of the free 
energy barrier. 

In conclusion, we have developed a modification of the single cluster Wolff algorithm for 
nematic systems which has enabled us to efficiently study the NI transition in Lebwohl- 
Lasher systems of sufficiently large size that finite size scaling is obeyed. As in the case of 
the cluster algorithms developed originally for ferromagnetic models, our algorithm allows us 
to overcome the critical slowing down associated with conventional single-flip Monte Carlo. 
The phase space can then be sampled efficiently near the transition as the system will flip 
readily between the ordered and disordered phases. We have also applied our algorithm to 
study the behavior of disclination loops in the transition region [[lj| and the efficiency of our 
algorithm should allow the study of many other interesting properties of the transition. 

We thank Prof. J. M. Kosterlitz for many helpful discussions and suggestions. This work 
was supported by the National Science Foundation under grant DMR-9873849. Computa- 
tional work in support of this research was performed at Brown University's Theoretical 
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FIGURES 

FIG. 1. (a) Illustration of the reflection operation -R(r), Eq. (|2|). The unit vector f is chosen 
randomly at the start of the algorithm. The reflection operation yields the new molecular orienta- 
tion a'i = R(r)<Ji as shown, (b) Illustration of the formation of a cluster of two molecules. Here we 
show the projections of the two molecular long axes ex, and Uj on the plane perpendicular to f , as 
well as the projections of the molecular long axes produced by the reflection operator R(r ) acting 
on these two molecules. Two molecules are likely to form a cluster if they each make an angle of 
approximately 45° with f and if the angle <f> between their projections is less than 90°. 

FIG. 2. Free energy, Eq. in units of e as a function of the energy per unit site E (also 
measured in units of e), for four different lattice sizes, L = 30 (•), 50 (A), 60 (o), 70 (*). The data 
for the three largest system sizes have been displaced vertically for the sake of clarity. 

FIG. 3. The free energy barrier height AF divided by L (measured in units of e over the lattice 
spacing) as a function of L (measured in units of the lattice spacing). The straight line fit for 
system sizes L > 35 indicates that AF obeys the finite size scaling relation, AF ~ L 2 , as expected 
for a first-order transition. 

FIG. 4. The transition temperature T C (L), (measured in units of e/kg) as a function of L~ 3 
(in units of 5 x 10 -5 cubic lattice spacings), for the eight system sizes shown in Fig. ||, showing 
the expected finite size scaling behavior (the straight line fit) given by Eq. (0) for system sizes 
L > 35. The extrapolation of this line to infinite system size yields an estimate of the transition 
temperature (indicated by the arrow) in the thermodynamic limit. 
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